To quantitatively examine the efficacy of vegetation restoration in drylands globally.
Study-level viz to document patterns in exclusions primarily and the relatie frequenices, at the study level, of major categories of evidence.
#study data####
library(tidyverse)
studies <- read_csv("data/studies.csv")
studies
## # A tibble: 278 x 18
## ID title technique data region exclude rationale observations
## <dbl> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 152 Shor… seeding,… expe… Africa no <NA> <NA>
## 2 180 Rest… chemical… App.… Africa no <NA> <NA>
## 3 229 Infl… soil see… fiel… Africa no <NA> <NA>
## 4 230 Acti… planting fiel… Africa no <NA> <NA>
## 5 255 The … grazing … fiel… Africa no <NA> <NA>
## 6 262 Reve… seeding,… eper… Africa no <NA> <NA>
## 7 263 The … phytogen… fiel… Africa no <NA> <NA>
## 8 264 Eval… seeding,… fiel… Africa no <NA> <NA>
## 9 271 Patc… natural … fiel… Africa no <NA> <NA>
## 10 4 Fact… natural … App.… Africa no <NA> <NA>
## # ... with 268 more rows, and 10 more variables: disturbance <chr>,
## # system <chr>, goal <chr>, intervention <chr>, paradigm <chr>,
## # grazing <chr>, hypothesis <chr>, soil <chr>, benchmark <chr>,
## # notes <chr>
#quick look at rationale needed
exclusions <- studies %>%
filter(exclude == "yes")
#quick look at studies with paradigms
evidence <- studies %>%
filter(exclude == "no")
#library(skimr)
#skim(evidence)
#study-level viz#####
#exclusions
ggplot(exclusions, aes(rationale, fill = region)) +
geom_bar() +
coord_flip() +
labs(x = "rational for exclusion", y = "frequency") +
scale_fill_brewer(palette = "Paired")
ggplot(evidence, aes(disturbance, fill = paradigm)) +
geom_bar(na.rm = TRUE) +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(y = "frequency")
ggplot(evidence, aes(region, fill = paradigm)) +
geom_bar(na.rm = TRUE) +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(y = "frequency")
ggplot(evidence, aes(data, fill = paradigm)) +
geom_bar(na.rm = TRUE) +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(y = "frequency")
ggplot(evidence, aes(system, fill = paradigm)) +
geom_bar(na.rm = TRUE) +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(y = "frequency")
ggplot(evidence, aes(goal, fill = paradigm)) +
geom_bar(na.rm = TRUE) +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(x = "outcome", y = "frequency")
#step 1 models####
#paradigm
derived.evidence <- evidence %>%
group_by(technique, data, region, disturbance, goal, paradigm) %>% summarise(n = n())
#active-passive split
m <- glm(n~paradigm, family = poisson, derived.evidence)
anova(m, test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 167 9.9147
## paradigm 1 0.045115 166 9.8696 0.8318
#region
m1 <- glm(n~paradigm*region, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m1, test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 167 9.9147
## paradigm 1 0.045115 166 9.8696 0.8318
## region 6 0.301367 160 9.5682 0.9995
## paradigm:region 6 0.213627 154 9.3546 0.9998
#outcome
m2 <- glm(n~paradigm*goal, family = poisson, derived.evidence)
#m1
#summary(m1)
anova(m2, test="Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: n
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 167 9.9147
## paradigm 1 0.045115 166 9.8696 0.8318
## goal 6 0.240941 160 9.6287 0.9997
## paradigm:goal 4 0.301480 156 9.3272 0.9897
#even split between active and passive evidence by all key categories
A summary of sort process using PRISMA.
library(PRISMAstatement)
prisma(found = 1504,
found_other = 5,
no_dupes = 1039,
screened = 1039,
screen_exclusions = 861,
full_text = 178,
full_text_exclusions = 100,
qualitative = 78,
quantitative = 78,
width = 800, height = 800)
Check data and calculate necessary measures.
data <- read_csv("data/data.csv")
data <- data %>%
mutate(lrr = log(mean.t/mean.c), rii = ((mean.t-mean.c)/(mean.t + mean.c)), var.es = ((sd.t^2/n.t*mean.t^2) + (sd.c^2/n.c*mean.c^2)))
#data
#consider adding some other effect size measures and/or study-level data too
Explore summary level data of all data. Explore aggregation levels that support the most reasonable data structure and minimize non-independence issues.
#evidence map####
require(maps)
world<-map_data("world")
map<-ggplot() + geom_polygon(data=world, fill="gray50", aes(x=long, y=lat, group=group))
map + geom_point(data=data, aes(x=long, y=lat)) +
labs(x = "longitude", y = "latitude") #render a literal map, i.e. evidence map, of where we study the niche in deserts globally
#add in levels and color code points on map####
map + geom_point(data=data, aes(x=long, y=lat, color = paradigm)) +
scale_color_brewer(palette = "Paired") +
labs(x = "longitude", y = "latitude", color = "")
#aggregation####
se <- function(x){
sd(x)/sqrt(length(x))
}
data.simple <- data %>%
group_by(study.ID, paradigm, technique, measure.success) %>%
summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))
main.data <- data %>%
group_by(study.ID, paradigm, intervention, outcome) %>%
summarise(n = n(), mean.lrr = mean(lrr), mean.rii = mean(rii), mean.var = mean(var.es))
#EDA data####
simple.data <- data %>% group_by(study.ID, paradigm, technique, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
simple.data <- na.omit(simple.data)
parad.data <- data %>% group_by(study.ID, paradigm) %>% summarise(mean.rii = mean(rii), error = se(rii))
parad.data <- na.omit(parad.data)
tech.data <- data %>% group_by(study.ID, technique) %>% summarise(mean.rii = mean(rii), error = se(rii))
tech.data <- na.omit(tech.data)
success.data <- data %>% group_by(study.ID, measure.success) %>% summarise(mean.rii = mean(rii), error = se(rii))
success.data <- na.omit(success.data)
#active
active <- data %>%
filter(paradigm == "active")
#viz for aggregation####
ggplot(na.omit(data.simple), aes(technique, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired")
ggplot(na.omit(data.simple), aes(measure.success, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired")
#better
ggplot(main.data, aes(intervention, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(fill = "")
ggplot(main.data, aes(outcome, n, fill = paradigm)) +
geom_bar(stat = "identity") +
coord_flip() +
scale_fill_brewer(palette = "Paired") +
labs(fill = "")
Exploratory data analyses to understand data and QA/QC using Rii.
#library(plotrix) #for quick s.e. calculations sometimes needed for data tidy step
#library(meta) #nice package for most meta-statistics
#assign model (typically a nice meta. function from one of several packages such as meta, metafor, or netmeta)
#need to write a quick function here to iterate and/or up split-map from broom
#EDA####
#rii####
#paradigm
#m <- metagen(mean.rii, error, studlab = study.ID, byvar = paradigm, data = simple.data) #fit generic meta-analysis to an object
#summary(m)
#no difference in random effects model by paradigm but fixed yes. Hetereogeneity is significantly different so a. fixed not representative and b. need a better model
#forest(m, xlim="symmetric", plotwidth=unit(1, "cm"))
#forest(m, sortvar = paradigm)
#radial(m)
#metabias(m, method = "linreg")
#technique
#m <- metagen(mean.rii, error, studlab = study.ID, byvar = technique, data = tech.data) #fit generic meta-analysis to an object
#m
#forest(m, sortvar = technique)
#radial(m)
#metabias(m, method = "linreg")
#measure
#m <- metagen(mean.rii, error, studlab = study.ID, byvar = measure.success, data = success.data) #fit generic meta-analysis to an object
#summary(m)
#forest(m, sortvar = measure.success)
#radial(m)
#metabias(m, method = "linreg")
#check data coding for passive intervention and outcome vectors too - not just active?
#explore above with intervention and outcome vectors
#do with ID as level not study.ID
#lrr
#rii anew with all data?
#m <- metagen(rii, var.es, studlab = ID, byvar = paradigm, data = data)
#summary(m)
#forest(m, sortvar = paradigm)
#radial(m)
#metabias(m, method = "linreg")
#m <- metagen(rii, var.es, studlab = study.ID, byvar = intervention, data = active)
#summary(m)
#forest(m, sortvar = paradigm)
#m <- metagen(rii, var.es, studlab = ID, byvar = outcome, data = na.exclude(data))
#summary(m)
#however repeat above for active only - for intervention and outcome - makes more sense
Meta and conventional statistical models to explore relative efficacy.
#effect sizes plots
ggplot(main.data, aes(paradigm, mean.rii, color = intervention)) +
geom_point(position = position_dodge(width = 0.5)) +
ylim(c(-1,1)) +
labs(x = "", y = "mean rii", color = "") +
coord_flip() +
geom_errorbar(aes(ymin=mean.rii-mean.var, ymax=mean.rii+mean.var), width=.05, position = position_dodge(width = 0.5)) +
geom_hline(yintercept = 0, colour="grey", linetype = "longdash")
ggplot(simple.data, aes(paradigm, mean.rii, color = technique)) +
geom_point(position = position_dodge(width = 0.5)) +
ylim(c(-1,1)) +
labs(x = "", y = "mean rii", color = "") +
coord_flip() +
geom_errorbar(aes(ymin=mean.rii-error, ymax=mean.rii+error), width=.05, position = position_dodge(width = 0.5)) +
geom_hline(yintercept = 0, colour="grey", linetype = "longdash")
ggplot(main.data, aes(paradigm, mean.rii, color = outcome)) +
geom_point(position = position_dodge(width = 0.5)) +
ylim(c(-1,1)) +
labs(x = "", y = "mean rii") +
coord_flip() +
geom_errorbar(aes(ymin=mean.rii-mean.var, ymax=mean.rii+mean.var), width=.05, position = position_dodge(width = 0.5)) +
geom_hline(yintercept = 0, colour="grey", linetype = "longdash")
#ggplot(simple.data, aes(paradigm, mean.rii, color = measure.success)) +
#geom_point(position = position_dodge(width = 0.5)) +
#ylim(c(-1,1)) +
#labs(x = "", y = "mean rii") +
#coord_flip() +
#geom_errorbar(aes(ymin=mean.rii-error, ymax=mean.rii+error), width=.05, position = position_dodge(width = 0.5)) +
#geom_hline(yintercept = 0, colour="grey", linetype = "longdash")
ggplot(active, aes(intervention, rii, color = outcome)) +
geom_point(position = position_dodge(width = 0.5)) +
ylim(c(-1,1)) +
labs(x = "", y = "rii") +
coord_flip() +
geom_errorbar(aes(ymin=rii-var.es, ymax=rii+var.es), width=.05, position = position_dodge(width = 0.5)) +
geom_hline(yintercept = 0, colour="grey", linetype = "longdash")
ggplot(active, aes(intervention, rii, color = outcome)) +
geom_boxplot() +
ylim(c(-1,1)) +
labs(x = "", y = "rii") +
coord_flip() +
geom_errorbar(aes(ymin=rii-var.es, ymax=rii+var.es), width=.05, position = position_dodge(width = 0.5)) +
geom_hline(yintercept = 0, colour="grey", linetype = "longdash")
#simple.active <- active %>%
#group_by(study.ID, outcome) %>%
#summarise(n = n(), mean.rii = mean(rii), mean.var = #mean(var.es))
#t-tests if different from 0
#t-tests if mean rii different from mu = 0
tmu <- function(x)
{t.test(x, mu = 0, paired = FALSE, var.equal=FALSE, conf.level = 0.95)
}
#data %>%
#split(.$paradigm) %>%
#purrr::map(~tmu(.$rii))
#active %>%
#split(.$intervention) %>%
#purrr::map(~tmu(.$rii))
#active %>%
#split(.$outcome) %>%
#purrr::map(~tmu(.$rii))
#meta-stats for active versus passive
#meta-stats for differences between interventions and outcomes for active test studies